setwd("I:/R Project/ConsensusClusterPlus")
###https://www.bioconductor.org/packages/release/data/experiment/html/ALL.html
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
library("BiocManager")
BiocManager::install("ALL")
BiocManager::install("ConsensusClusterPlus")


#TCGA数据
setwd("D:/R/table")
load("dataTCGANorm.rda")
#我做的时候只做了前10000个，后面这个方括号丢掉就好
TCGA <- dataTCGANorm[1:10000, 1:670]
rm(dataTCGANorm)
class(TCGA)
#稍微处理下数据，不用改
mads=apply(TCGA,1,mad)
TCGA=TCGA[rev(order(mads))[1:5000],]
# agglomerative hierarchical clustering algorithm using Pearson correlation distance归一化操作
TCGA = sweep(TCGA,1, apply(TCGA,1,median,na.rm=T))
#一键出图，具体后台做了个啥表参照师兄发我的手册ConsensusClusterPlus (Tutorial)中
#3.3 Generating cluster and item consensus部分
library(ConsensusClusterPlus)
results = ConsensusClusterPlus(TCGA,#分析矩阵
                               maxK=6,#最大聚类的数目，建议20
                               reps=1000,#重抽样的次数，建议1000
                               pItem=0.8,#sample重抽样的比例
                               pFeature=1,#gene重抽样的比例
                               title="D:/neurosurgery/R Project/ConsensusClusterPlus/plots",#保存路径
                               clusterAlg="hc",#聚类算法，此处为agglomerative hierarchical clustering algorithm
                               distance="pearson",#计算距离的方法，此处为pearson相关系数
                               seed=1262118388.71279,#设置随机种子，方便重复
                               plot="pdf"#结果图片的导出类型
)

icl <- calcICL(results,title ="D:/neurosurgery/R Project/ConsensusClusterPlus/xiangya",plot = 'pdf')

##下面基本是郑大一附院共识聚类的代码，代码链接https://github.com/Zaoqu-Liu/IRLS/blob/main/Consensus%20cluster.R
#PAC = Proportion of ambiguous clustering 模糊聚类比例，用于进一步确定最佳聚类数
Kvec = 2:9
x1 = 0.1; x2 = 0.9 
PAC = rep(NA,length(Kvec)) 
names(PAC) = paste("K=",Kvec,sep="") 
for(i in Kvec){
  M = results[[i]]$consensusMatrix
  Fn = ecdf(M[lower.tri(M)])
  PAC[i-1] = Fn(x2) - Fn(x1)
}
optK = Kvec[which.min(PAC)]
optK

PAC <- as.data.frame(PAC)
PAC$K <- 2:9
library(ggplot2)
ggplot(PAC,aes(factor(K),PAC,group=1))+
  geom_line()+
  theme_bw(base_rect_size = 1.5)+
  geom_point(size=4,shape=21,color='darkred',fill='orange')+
  ggtitle('Proportion of ambiguous clustering')+
  xlab('Cluster number K')+ylab(NULL)+
  theme(axis.text = element_text(size=12),
        plot.title = element_text(hjust=0.5),
        axis.title = element_text(size=13))
ggsave(filename = 'D:/neurosurgery/R Project/ConsensusClusterPlus/plots/PAC.pdf',width = 3.8,height = 4)

## 保存分型信息,为每个样本加一个feature叫做cluster，如此即可和别的模块对接
clusterNum=2      
cluster=results[[clusterNum]][["consensusClass"]]
#sub即可和别的模块对接
sub <- data.frame(Sample=names(cluster),Cluster=cluster)
sub$Cluster <- paste0('C',sub$Cluster)
table(sub$Cluster)

head(sub)

#矩阵图改成分型图，更专业一点
my <- results[[2]][["ml"]]
library(pheatmap)
rownames(my) <- sub$Sample
colnames(my) <- sub$Sample
pheatmap(1-my,show_colnames = F,show_rownames = F,
         treeheight_row = 20,treeheight_col = 20,
         clustering_method = 'complete',
         color = colorRampPalette(c("white","#C75D30"))(50), 
         annotation_names_row = F,annotation_names_col = F,
         annotation_row = data.frame(Cluster=sub$Cluster,row.names = sub$Sample),
         annotation_col = data.frame(Cluster=sub$Cluster,row.names = sub$Sample),
         annotation_colors = list(Cluster=c('C2'='#B5739D','C1'='#4E8279')))
library(export)
graph2pdf(file='D:/neurosurgery/R Project/ConsensusClusterPlus/plots/pheatmap.pdf',width=5.5,height=4.5)

